module dbessels
implicit none
contains
function vdlogk(vnu,x,tol)
!Derivative of the logarithm of the modified Bessel function of the third (logK(vnu,x)) kind with respect to order (vnu)
!(see appendix C in Mencia and Sentana 2004)
use bessels
double precision vdlogk,vnu,x,tol,nu,errorent,res,f0,f2,parte1,parte2,nui
double precision, parameter:: pi=3.141592653589793108624468950438d0
    nu=dabs(vnu)
    errorent=nu-round(nu)
    if (nu > 10.0d0) then
        res=dlogkxl(nu,x)
    elseif (dabs(x).le. 12.0d0) then
        if (nu==0.0d0) then
            res=0.0d0
        elseif (errorent==0.0d0) then
            nui=round(nu)
            f0=besseli(-nui,x)-besseli(nui,x)
            f2=d2idnu2(-nui,x,tol/10.0d0)-d2idnu2(nui,x,tol/10.0d0)
            res=((f2+pi**2.0d0*f0)/(4.0d0*dcos(nui*pi)))/besselk(nu,x)
        elseif (abs(errorent)<.03d0) then
            nui=round(nu)
            f0= 0.0d0
            f2=(d2idnu2(-nui,x,tol/10.0d0)-d2idnu2(nui,x,tol/10.0d0))/(4.0d0*dcos(nui*pi))
			
			res=f0+f2&
			&+d2kdnu2(nui,x,tol/100.0d0)*(nu-nui)&
			&+.5*d3kdnu3(nui,x,tol/100.0d0)*((nu-nui)**2.0d0)&
			&+(1/6.0d0)*d4kdnu4(nui,x,tol/100.0d0)*((nu-nui)**3.0d0)
            res=res/besselk(nu,x)

        else
            parte1=(-didnu(-nu,x,tol)-didnu(nu,x,tol))*csc(pi*nu)*.5d0*pi
            parte2=pi*dcotan(pi*nu)*besselk(nu,x)
            res=(parte1-parte2)/besselk(nu,x)
        end if
	elseif (dabs(x)>100.0d0) then
	    res=dlogkdnularge(nu,x)
    else
        res=dkdnularge(nu,x)/besselk(nu,x)
    end if
if (vnu .ge. 0) then
    vdlogk=res
else
    vdlogk=-res
end if
end function vdlogk
!%%%%%%%%%%%%%%%%%%%%%% Auxiliary Programs %%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dlogkxl(nu,z)
!Derivative of logK(nu,z) with respect to order
!Valid when both nu and z are large
double precision nu,z,dlogkxl,coc,raiz,t,tp,u1,u1p,u2,u2p,res
coc=z/nu
raiz=dsqrt(1.0d0+(coc*coc))
t=1.0d0/raiz
tp=coc*coc/(nu*(raiz**3.0d0))

u1=(3.0d0*t-5.0d0*(t**3.0d0))/24.0d0
u1p=(3.0d0-15.0d0*(t**2.0d0))/24.0d0

u2=(81.0d0*(t**2.0d0)-462.0d0*(t**4.0d0)+385*(t**6.0d0))/1152.0d0
u2p=(162.0d0*t-1848.0d0*(t**3.0d0)+2310*(t**5.0d0))/1152.0d0

res=(u1-2.0d0*(u2/nu)+(-u1p*nu+u2p)*tp)/((nu**2.0d0)-u1*nu+u2)
dlogkxl=1.0d0+dlog(nu)-(.5d0/nu)-dlog(z)-raiz&
&+((coc**2.0d0)/raiz)+((coc**2.0d0)/(2.0d0*nu*(raiz**2.0d0)))&
&+dlog(1.0d0+raiz)-((coc**2.0d0)/(1.0d0+(coc**2.0d0)+raiz))+res

end function dlogkxl

function dkdnularge(nu,x)
!Derivative of K(nu,x) with respect to nu for x large
double precision nu,x,dkdnularge,coef,muu,sum1,sum2,sum3,sum4,sum5,sum6,dgamma
double precision, parameter:: pi=3.141592653589793108624468950438d0
coef=dsqrt(pi/(2.0d0*x))*dexp(-x)
muu=4.0d0*nu**2.0d0
sum1=nu/x
sum2=nu*(2.0d0*muu-10.0d0)/(16.0d0*x**2)
sum3=nu*(3.0d0*muu**2.0d0-70.0d0*muu+259.0d0)/(384.0d0*x**3)

sum4=8.0d0*nu*((muu-9.0d0)*(muu-25.0d0)*(muu-49.0d0)&
&+(muu-1.0d0)*(muu-25.0d0)*(muu-49.0d0)&
&+(muu-1.0d0)*(muu-9.0d0)*(muu-49.0d0)+(muu-1.0d0)*(muu-9.0d0)*(muu-25.0d0))/(24.0d0*(8.0d0*x)**4.0d0)

sum5=8.0d0*nu*((muu-9.0d0)*(muu-25.0d0)*(muu-49.0d0)*(muu-81.0d0)&
&+(muu-1.0d0)*(muu-25.0d0)*(muu-49.0d0)*(muu-81.0d0)+(muu-1.0d0)*(muu-9.0d0)*(muu-49.0d0)*(muu-81.0d0)&
&+(muu-1.0d0)*(muu-9.0d0)*(muu-25.0d0)*(muu-81.0d0)&
&+(muu-1.0d0)*(muu-9.0d0)*(muu-25.0d0)*(muu-49.0d0))/(dgamma(6.0d0)*(8.0d0*x)**5.0d0)


sum6=8.0d0*nu*(6.0d0*muu**5.0d0-5.0d0*286.0d0*muu**4.0d0+4.0d0*28743.0d0*muu**3.0d0&
&-3.0d0*1234948.0d0*muu**2.0d0+2.0d0*21967231.0d0*muu-128816766.0d0)&
&/(dgamma(7.0d0)*(8.0d0*x)**6.0d0)
dkdnularge=coef*(sum1+sum2+sum3+sum4+sum5+sum6)
end function dkdnularge


function dlogkdnularge(nu,x)
!Derivative of logK(nu,x) with respect to nu for x large
double precision nu,x,dlogkdnularge,p1,p2,p3,p4,p5,m,den,ix,num,dp1,dp2,dp3,dp4,dp5
double precision, parameter:: pi=3.141592653589793108624468950438d0

m=4.0d0*nu*nu
ix=1.0d0/(8.0d0*x)
p1=m-1.0d0
p2=p1*(m-9.0d0)
p3=p2*(m-25.0d0)
p4=p3*(m-49.0d0)
p5=p4*(m-81.0d0)

dp1=8.0d0*nu
dp2=dp1*(m-9.0d0)+p1*8.0d0*nu
dp3=dp2*(m-25.0d0)+p2*8.0d0*nu
dp4=dp3*(m-49.0d0)+p3*8.0d0*nu
dp5=dp4*(m-81.0d0)+p4*8.0d0*nu

num=dp1*ix+(dp2*(ix**2.0d0)/2.0d0)+(dp3*(ix**3.0d0)/6.0d0)+(dp4*(ix**4.0d0)/24.0d0)+(dp5*(ix**5.0d0)/120.0d0)
den=1+p1*ix+(p2*(ix**2.0d0)/2.0d0)+(p3*(ix**3.0d0)/6.0d0)+(p4*(ix**4.0d0)/24.0d0)+(p5*(ix**5.0d0)/120.0d0)

dlogkdnularge=num/den
end function dlogkdnularge

!%%%%%

function dkdnu(nu,x,tol)
!Derivative of  K(nu,x) with respect to nu
use bessels
double precision nu,x,tol,dkdnu,errorent,nui,f0,f2,parte1,parte2,res
double precision, parameter:: pi=3.141592653589793108624468950438d0
errorent=nu-round(nu)
if (nu > 10.0d0) then
    res=dlogkxl(nu,x)*besselk(nu,x)
elseif (abs(x).le. 12.0d0) then
    if (nu==0.0d0) then
        res=0
    elseif (errorent==0) then
        nui=round(nu)
        f0=besseli(-nui,x)-besseli(nui,x)
        f2=d2idnu2(-nui,x,tol/10.0d0)-d2idnu2(nui,x,tol/10.0d0)
        res=((f2+pi**2.0d0*f0)/(4.0d0*dcos(nui*pi)))
    
    elseif (abs(errorent)<.1d0) then
        nui=round(nu)
        f0=besseli(-nui,x)-besseli(nui,x)
        f2=d2idnu2(-nui,x,tol/10.0d0)-d2idnu2(nui,x,tol/10.0d0)
        res=((f2+pi**2.0d0*f0)/(4.0d0*dcos(nui*pi)))&
		&+d2kdnu2(nui,x,tol/10.0d0)*(nu-nui)+.5d0*d3kdnu3(nui,x,tol/10.0d0)*((nu-nui)**2.0d0)&
		&+(1.0d0/6.0d0)*d4kdnu4(nui,x,tol/10.0d0)*((nu-nui)**3.0d0)
    else
        parte1=(-didnu(-nu,x,tol)-didnu(nu,x,tol))*csc(pi*nu)*.5d0*pi
        parte2=pi*dcotan(pi*nu)*besselk(nu,x)
        res=parte1-parte2
    end if
else
    res=dkdnularge(nu,x)
end if
dkdnu=res
end function dkdnu
!%%%%%

function d2kdnu2(nu,x,tol)
!Second derivative of K(nu,x) with respect to nu
use bessels
double precision d2kdnu2,nu,x,tol,zer,nui,arg1,arg2,res
double precision, parameter:: pi=3.141592653589793108624468950438d0


zer=round(nu)-nu
zer=abs(zer)

if (zer==0) then
    nui=round(nu)
    arg1=(-d3idnu3(-nui,x,tol/4.0d0)-d3idnu3(nui,x,tol/4.0d0))/(6.0d0*dcos(nui*pi))
    arg2=(pi**2.0d0)*(-didnu(-nui,x,tol/4.0d0)-didnu(nui,x,tol/4.0d0))/(3.0d0*dcos(nui*pi))
    d2kdnu2=arg1+arg2-((pi**2.0d0)*besselk(nui,x)/3.0d0)    
elseif (zer < .1d0) then
    nui=round(nu)
    arg1=(-d3idnu3(-nui,x,tol/4.0d0)-d3idnu3(nui,x,tol/4.0d0))/(6.0d0*dcos(nui*pi))
    arg2=(pi**2.0d0)*(-didnu(-nui,x,tol/4.0d0)-didnu(nui,x,tol/4.0d0))/(3.0d0*dcos(nui*pi))
    res=arg1+arg2-((pi**2.0d0)*besselk(nui,x)/3.0d0)
    d2kdnu2=res+d3kdnu3(nui,x,tol)*(nu-nui)+.5d0*d4kdnu4(nui,x,tol)*((nu-nui)**2.0d0)
else
    d2kdnu2=(d2idnu2(-nu,x,tol/3.0d0)-d2idnu2(nu,x,tol/3.0d0)&
	&-4.0d0*dcos(nu*pi)*dkdnu(nu,x,tol/3.0d0)&
	&+(pi**2.0d0)*(besseli(-nu,x)-besseli(nu,x)))*pi/(2.0d0*dsin(nu*pi))
end if

end function d2kdnu2
!%%%%%

function d3kdnu3(nu,x,tol)
!Third derivative of K(nu,x) with respect to nu
use bessels
double precision d3kdnu3,nu,x,tol,zer,nui,f0,f1,f2,f4,res,arg1,arg2
double precision, parameter:: pi=3.141592653589793108624468950438d0

zer=round(nu)-nu
zer=dabs(zer)

if (zer==0.0d0) then
    nui=round(nu)
    f0=besseli(-nui,x)-besseli(nui,x)
    f1=-didnu(-nui,x,tol)-didnu(nui,x,tol)
    f2=d2idnu2(-nui,x,tol)-d2idnu2(nui,x,tol)
    f4=d4idnu4(-nui,x,tol)-d4idnu4(nui,x,tol)
    d3kdnu3=(f4+2.0d0*pi**2.0d0*f2+14.0d0*pi**3.0d0*dtan(nui*pi)*f1-6.0d0*pi**4.0d0*f0)&
	&/(8.0d0*dcos(nui*pi))    
elseif (zer < .07d0) then
    nui=round(nu)
    f0=besseli(-nui,x)-besseli(nui,x)
    f1=-didnu(-nui,x,tol)-didnu(nui,x,tol)
    f2=d2idnu2(-nui,x,tol)-d2idnu2(nui,x,tol)
    f4=d4idnu4(-nui,x,tol)-d4idnu4(nui,x,tol)
    res=(f4+2.0d0*pi**2.0d0*f2+14.0d0*pi**3.0d0*dtan(nui*pi)*f1-6.0d0*pi**4.0d0*f0)/(8.0d0*dcos(nui*pi))
    d3kdnu3=res+d4kdnu4(nui,x,tol/10.0d0)*(nu-nui)
else
    arg1=.5d0*pi*(-d3idnu3(-nu,x,tol/4.0d0)-d3idnu3(nu,x,tol/4.0d0))/dsin(pi*nu)
    arg2=-3.0d0*pi*dcotan(pi*nu)*d2kdnu2(nu,x,tol/4.0d0)+3.0d0*(pi**2.0d0)*dkdnu(nu,x,tol/4.0d0)&
	&+(pi**3.0d0)*dcotan(pi*nu)*besselk(nu,x)
    
    d3kdnu3=arg1+arg2    
end if
end function d3kdnu3
!%%%%%

function d4kdnu4(nu,x,tol)
!Fourth derivative of K(nu,x) with respect to nu
use bessels
double precision d4kdnu4,nu,x,tol,zer,nui,f0,f1,f2,f3,f5,arg1,arg2
double precision, parameter:: pi=3.141592653589793108624468950438d0


zer=round(nu)-nu
zer=abs(zer)

if (zer==0.0d0) then
    nui=round(nu)
    f0=besseli(-nui,x)-besseli(nui,x)
    f1=-didnu(-nui,x,tol/10.0d0)-didnu(nui,x,tol/10.0d0)
    f2=d2idnu2(-nui,x,tol/10.0d0)-d2idnu2(nui,x,tol/10.0d0)
    f3=-d3idnu3(-nui,x,tol/10.0d0)-d3idnu3(nui,x,tol/10.0d0)
    f5=-d5idnu5(-nui,x,tol/10.0d0)-d5idnu5(nui,x,tol/10.0d0)
    d4kdnu4=((1.5d0*f5-10.0d0*pi**2.0d0*f3-4.0d0*pi**4.0d0*f1)/(15.0d0*dcos(nui*pi)))&
	&+6.0d0*pi**2.0d0*d2kdnu2(nui,x,tol/10.0d0)-pi**4.0d0*besselk(nui,x)    
    
elseif (zer < .08d0) then
    nui=round(nu)
    f0=besseli(-nui,x)-besseli(nui,x)
    f1=-didnu(-nui,x,tol/10.0d0)-didnu(nui,x,tol/10.0d0)
    f2=d2idnu2(-nui,x,tol/10.0d0)-d2idnu2(nui,x,tol/10.0d0)
    f3=-d3idnu3(-nui,x,tol/10.0d0)-d3idnu3(nui,x,tol/10.0d0)
    f5=-d5idnu5(-nui,x,tol/10.0d0)-d5idnu5(nui,x,tol/10.0d0)
    d4kdnu4=((1.5d0*f5-10.0d0*pi**2.0d0*f3-4.0d0*pi**4.0d0*f1)/(15.0d0*dcos(nui*pi)))&
	&+6.0d0*pi**2.0d0*d2kdnu2(nui,x,tol/10.0d0)-pi**4.0d0*besselk(nui,x)    
    
else
    arg1=.5*pi*(d4idnu4(-nu,x,tol/10.0d0)-d4idnu4(nu,x,tol/10.0d0))/dsin(pi*nu)
	arg2=-4.0d0*pi*dcotan(pi*nu)*d3kdnu3(nu,x,tol/10.0d0)&
	&+6.0d0*(pi**2.0d0)*d2kdnu2(nu,x,tol/10.0d0)&
	&+4.0d0*(pi**3.0d0)*dcotan(pi*nu)*dkdnu(nu,x,tol/10.0d0)-pi**4.0d0*besselk(nu,x)
    d4kdnu4=arg1+arg2    
end if
end function d4kdnu4
!%%%%%

function didnu(nu,x,tol)
!Derivative of the Bessel function of the third kind I(nu,x) with respect to order(nu) (nu is integer)
use bessels
double precision didnu,nu,x,tol,res1,res2,contk,obj,orden,iorden,coc,fac,ad,dgamma,dpsi
double precision, parameter:: pi=3.141592653589793108624468950438d0
res1=besseli(nu,x)*dlog(.5*x)
res2=0.0d0
contk=0.0d0
obj=1.0d5+tol
do while (obj>tol)
    orden=nu+contk+1.0d0
    if (orden > 0.0d0) then
        coc=dpsi(orden)/dgamma(orden)
      else
        iorden=1.0d0-orden
        coc=(dpsi(iorden)*dsin(pi*orden)-pi*dcos(pi*orden))*(dgamma(iorden)/pi)
    end if
    fac=(.25d0*(x**2.0d0))**contk
    ad=coc*fac/dgamma(contk+1.0d0)
    obj=ad*((.5d0*x)**nu)
    res2=res2+obj
    contk=contk+1.0d0
    obj=dabs(obj)
end do

didnu=res1-res2
end function didnu
!%%%%%

function d2idnu2(nu,x,tol)
!Second derivative of I(nu,x) with respect to nu (nu is integer)
use bessels
use polygam
double precision d2idnu2,nu,x,tol,a,res1,res2,contk,orden,coc,iorden,obj,dgamma,dpsi
double precision, parameter:: pi=3.141592653589793108624468950438d0

a=dlog(.5d0*x)
res1=-besseli(nu,x)*(a**2)+2*a*didnu(nu,x,tol/10.0d0)
res2=0.0d0
contk=0.0d0
obj=1.0d5
do while (obj>(tol/2.0d0)) 
    orden=nu+contk+1
    if (orden > 0.0d0) then
        coc=(ppsi(1.0d0,orden)-((dpsi(orden))**2.0d0))/dgamma(orden)
      else
        iorden=1-orden
        coc=(1.0d0/pi)*dgamma(iorden)*(((pi**2)-ppsi(1.0d0,iorden)&
		&-(dpsi(iorden)**2.0d0))*dsin(pi*orden)+2.0d0*pi*dpsi(iorden)*dcos(pi*orden))
        
     end if
    obj=((.5d0*x)**nu)*coc*((.25d0*(x**2.0d0))**contk)/dgamma(contk+1.0d0)
    res2=res2+obj
    contk=contk+1.0d0
    obj=dabs(obj)
end do


d2idnu2=res1-res2
end function d2idnu2
!%%%%%
function d3idnu3(nu,x,tol)
!Third derivative of I(nu,x) with respect to nu (nu is integer)
use bessels
use polygam
double precision d3idnu3, nu,x, tol,a,res1,res2,contk,obj,orden,iorden,arg1,arg2,coc,dpsi,dgamma
double precision, parameter:: pi=3.141592653589793108624468950438d0
a=dlog(.5d0*x)
res1=3.0d0*a*d2idnu2(nu,x,tol/10.0d0)-3.0d0*(a**2.0d0)*didnu(nu,x,tol/10.0d0)+(a**3.0d0)*besseli(nu,x)
res2=0.0d0

contk=0.0d0
obj=1.0d5
do while (obj>(tol/4.0d0))
    orden=nu+contk+1.0d0
    if (orden > 0.0d0) then
        coc=((dpsi(orden)**3.0d0)-3.0d0*ppsi(1.0d0,orden)*dpsi(orden)+ppsi(2.0d0,orden))/dgamma(orden)
      else
        iorden=1.0d0-orden
        arg1=(dpsi(iorden)**3.0d0)-3.0d0*(pi**2.0d0)*dpsi(iorden)&
		&+3.0d0*dpsi(iorden)*ppsi(1.0d0,iorden)+ppsi(2.0d0,iorden)
        arg2=-3.0d0*pi*dcos(pi*orden)*(dpsi(iorden)**2.0d0)&
		&+(pi**3.0d0)*dcos(pi*orden)-3.0d0*pi*dcos(pi*orden)*ppsi(1.0d0,iorden)
        coc=(1.0d0/pi)*dgamma(iorden)*(arg1*dsin(pi*orden)+arg2)
    end if
    obj=((.5d0*x)**nu)*coc*((.25d0*(x**2.0d0))**contk)/dgamma(contk+1.0d0)
    res2=res2+obj
    contk=contk+1.0d0
    obj=dabs(obj)
end do
 

d3idnu3=res1-res2
end function d3idnu3
!%%%%%

function d4idnu4(nu,x,tol)
!Fourth derivative of I(nu,x) with respect to nu (nu is integer)
use bessels
use polygam
double precision d4idnu4,nu,x,tol,a,res1,res2,contk,obj,orden,iorden,coc,p0,p1,arg1,arg2,dpsi,dgamma
double precision, parameter:: pi=3.141592653589793108624468950438d0
a=dlog(.5d0*x)
res1=4.0d0*a*d3idnu3(nu,x,tol/10.0d0)-6.0d0*(a**2.0d0)*d2idnu2(nu,x,tol/10.0d0)&
&+4.0d0*(a**3.0d0)*didnu(nu,x,tol/10.0d0)-(a**4.0d0)*besseli(nu,x)
res2=0.0d0

contk=0
obj=1.0d5
do while (obj>(tol/4.0d0))
    orden=nu+contk+1.0d0
    if (orden > 0.0d0) then
        coc=((-dpsi(orden)**4.0d0)+6.0d0*(dpsi(orden)**2.0d0)*ppsi(1.0d0,orden)&
		&-4.0d0*dpsi(orden)*ppsi(2.0d0,orden)-3.0d0*(ppsi(1.0d0,orden)**2.0d0)&
		&+ppsi(3.0d0,orden))/dgamma(orden)
      else
        iorden=1.0d0-orden
        p0=dpsi(iorden)
        p1=ppsi(1.0d0,iorden)
        
        arg1=(1.0d0/pi)*dgamma(iorden)*(-(p0**4.0d0)+6.0d0*pi**2.0d0*p0**2.0d0-6.0d0*p0**2*p1&
		&-4.0d0*p0*ppsi(2.0d0,iorden)-3.0d0*p1**2.0d0-ppsi(3.0d0,iorden)+6.0d0*pi**2.0d0*p1-pi**4.0d0)
        arg2=dgamma(iorden)*(4.0d0*p0**3.0d0-4.0d0*pi**2.0d0*p0+12.0d0*p0*p1+4.0d0*ppsi(2.0d0,iorden))
        coc=arg1*dsin(pi*orden)+arg2*dcos(pi*orden)
    end if
    obj=((.5d0*x)**nu)*coc*((.25d0*(x**2.0d0))**contk)/dgamma(contk+1.0d0)
    res2=res2+obj
    contk=contk+1.0d0
    obj=dabs(obj)
end do
d4idnu4=res1-res2
end function d4idnu4
!%%%%

function d5idnu5(nu,x,tol)
!Fifth derivative of I(nu,x) with respect to nu (nu is integer)
use bessels
use polygam
double precision d5idnu5,nu,x,tol,a,res1,res2,contk,obj,orden,coc,iorden,p0,p1,p2,arg11&
&,arg12,arg13,arg1,arg21,arg2,dgamma,dpsi
double precision, parameter:: pi=3.141592653589793108624468950438d0
a=dlog(.5d0*x)
res1=5.0d0*a*d4idnu4(nu,x,tol/10.0d0)&
&-10.0d0*(a**2.0d0)*d3idnu3(nu,x,tol/10.0d0)+10.0d0*(a**3.0d0)*d2idnu2(nu,x,tol/10.0d0)&
&-5.0d0*(a**4.0d0)*didnu(nu,x,tol/10.0d0)+(a**5.0d0)*besseli(nu,x)
res2=0.0d0

contk=0.0d0
obj=1.0d5
do while (obj>(tol/4.0d0))
    orden=nu+contk+1.0d0
    if (orden > 0.0d0) then
        coc=((dpsi(orden)**5.0d0)-10.0d0*(dpsi(orden)**3.0d0)*ppsi(1.0d0,orden)&
		&+10.0d0*(dpsi(orden)**2.0d0)*ppsi(2.0d0,orden)+15.0d0*dpsi(orden)*(ppsi(1.0d0,orden)**2.0d0)&
		&-5.0d0*dpsi(orden)*ppsi(3.0d0,orden)&
		&-10.0d0*ppsi(1.0d0,orden)*ppsi(2.0d0,orden)+ppsi(4.0d0,orden))/dgamma(orden)
      else
        iorden=1.0d0-orden
        p0=dpsi(iorden)
        p1=ppsi(1.0d0,iorden)
        p2=ppsi(2.0d0,iorden)
        arg11=p0**5.0d0-10.0d0*pi**2.0d0*p0**3.0d0+10.0d0*p0**3.0d0*p1+10.0d0*p0**2.0d0*p2
        arg12=15.0d0*p0*p1**2.0d0+5.0d0*p0*ppsi(3.0d0,iorden)+5.0d0*pi**4.0d0*p0
        arg13=-30.0d0*pi**2.0d0*p0*p1+10.0d0*p1*p2-10.0d0*pi**2.0d0*p2+ppsi(4.0d0,iorden)
        arg1=(1.0d0/pi)*dgamma(iorden)*(arg11+arg12+arg13)*dsin(pi*orden)
        
        arg21=-5.0d0*p0**4.0d0+10.0d0*pi**2.0d0*p0**2.0d0-30.0d0*p0**2.0d0*p1&
		&-20.0d0*p0*p2-15.0d0*p1**2.0d0+10.0d0*pi**2.0d0*p1-5.0d0*ppsi(3.0d0,iorden)-pi**4.0d0
        arg2=dgamma(iorden)*arg21*dcos(pi*orden)
        coc=arg1+arg2
    end if
    obj=((.5d0*x)**nu)*coc*((.25d0*(x**2.0d0))**contk)/dgamma(contk+1.0d0)
    res2=res2+obj
    contk=contk+1.0d0
    obj=dabs(obj)
end do


d5idnu5=res1-res2
end function d5idnu5
!%%%%

function round(x)
double precision x,round
if ((x-floor(x)).le.(ceiling(x)-x)) then
	round=dble(floor(x))
else
	round=dble(ceiling(x))
end if
end function round

function csc(x)
double precision csc,x
csc=1.0d0/dsin(x)
end function csc
end module dbessels